基因组信息学实验

发表于 | 更新于

本学期总共六个实验，先后完成了对酵母YJM320的全基因组的模拟测序、序列组装、同源搜索、从头预测、启动子（TATA box）预测、全基因组数据可视化。

摘要

本学期总共六个实验，先后完成了对酵母YJM320的全基因组的模拟测序、序列组装、同源搜索、从头预测、启动子（TATA box）预测、全基因组数据可视化。
具体内容如下：利用art illumina工具对原始的酵母YJM320的全基因组进行模拟建模，利用fastqc对建模结果进行质控分析（质控结果很好）。然后对建模结果利用SOAPdenovo进行序列组装，组装结果为Scaffold265个，其N50 258107bp；contig大于100bp的共43082个，其N50为273bp。结果利用blastn和quast同原始基因组进行比对，两种方法计算覆盖度分别为75.94%和87.6%，利用IGV工具观察基因组，发现在12号染色体中有一段rRNA重复序列无法组装，这也是影响最终覆盖率的主要原因。
分别采用同源搜索（利用tblastn比对）和从头测序（Augustus）的方法构建全基因组，对所得结果同原基因组进行gffcompare比对，结果为从头测序结果优于同源建模。IGV是数据可视化的工具，利用它可以查看同源搜索和从头预测的结果。
采用HMM结合bootstrap预测全基因组中的TATA box分布，绘制打分图和ROC曲线图，有70%匹配的TATA box下游5000bp内有发现基因。

一、材料和方法

1.1、分析平台

1.1.1 硬件平台

(1).CPU: Intel(R) Core(TM) i7-7700HQ CPU @ 2.80GHz

(2). 内存：8 GB 1600 MHz DDR3

(3). 硬盘：磁盘 0 (C:) TOSHIBA THNSNK128GVN8 M.2 2280 128GB

磁盘 1 (D: E: F:) TOSHIBA MQ01ABD100 932 GB

1.1.2 操作系统

Windows10 、Ubuntu 虚拟机

1.1.3 分析软件

（含版本号）

编程工具：

R 3.5.1 Perl 5.26.1

分析工具：

ART 2.5.8 sratoolkit 2.9.4 fastqc SOAPdenovo

Quast 本地blast gffcompare Augustus

1.1.4 数据库资源

NCBI Genome：https://www.ncbi.nlm.nih.gov/genome/15?genome_assembly_id=341003

NCBI PubMed：https://www.ncbi.nlm.nih.gov/pubmed/

NCBI SRA：https://www.ncbi.nlm.nih.gov/sra/SRX257750[accn]]

UniProt：https://www.uniprot.org

EPD：https://epd.epfl.ch//index.php

1.2 研究对象

物种名称 Saccharomyces cerevisiae YJM320

Genbank Accession Number:GCA_000975885.2

1.3 方法

1.3.0 实验流程图

图 1 实验流程图

1.3.1 基因组测序模拟

分别下载YJM320的sra测序数据、gff注释数据和fna基因组数据，其中sra测序地址、gff数据地址，fna数据地址。使用art illumina软件进行序列模拟，模拟参数见表1，计算理论覆盖率。采用SRA Toolkit工具验证sra数据完整性，并从中提取fastq基因序列数据。

表 1 art illumina模拟参数设置

l（reads的长度） f（覆盖度） m（片段平均长度）
100 1 150
100 2 150
100 3 150
100 5 200
80 10 150

运行代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#art illumina序列模拟

art\_illumina -ss HS20 -sam -i GCA\_000975885.2\_Sc\_YJM320\_v1\_genomic.fna -p -l 100 -f 10 -m 150 -s 50 -o data/paired\_dat\_10\_

art\_illumina -ss HS20 -sam -i GCA\_000975885.2\_Sc\_YJM320\_v1\_genomic.fna -p -l 100 -f 5 -m 150 -s 10 -o data/paired\_dat\_5\_

art\_illumina -ss HS20 -sam -i GCA\_000975885.2\_Sc\_YJM320\_v1\_genomic.fna -p -l 100 -f 3 -m 150 -s 10 -o data/paired\_dat\_3\_

art\_illumina -ss HS20 -sam -i GCA\_000975885.2\_Sc\_YJM320\_v1\_genomic.fna -p -l 100 -f 2 -m 150 -s 10 -o data/paired\_dat\_2\_

art\_illumina -ss HS20 -sam -i GCA\_000975885.2\_Sc\_YJM320\_v1\_genomic.fna -p -l 100 -f 1 -m 150 -s 10 -o data/paired\_dat\_1\_

#SRA Toolkit数据校验

vdb-validate SRR800771.sra

#提取fastq数据

fastq-dump --split-3 ERR3139854.sra

1.3.2 序列组装

采用art illumina工具分别模拟180x和3000x两套数据，对数据进行fastqc质控分析，确认数据符合质控标准后，选择一套参考基因组数据，这里选用参考基因组数据。利用bowtie2分别进行两套数据同参考基因组的比对和samtools分析，这里我特别分析了测序深度方面的数据。利用SOAPdenovo软件对序列片段进行序列组装。利用网络平台QUAST（http://quast.bioinf.spbau.ru/）进行可视化及实际覆盖率计算。

运行代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
#art illumina工具分别模拟180x和3000x两套数据

./art\_illumina -ss HS20 -sam -i GCA\_000975675.2\_Sc\_YJM1083\_v1\_genomic.fna -p -l 100 -f 45 -m 3000 -s 20 -o paired\_dat\_m3000\_

./art\_illumina -ss HS20 -sam -i GCA\_000975675.2\_Sc\_YJM1083\_v1\_genomic.fna -p -l 100 -f 45 -m 180 -s 20 -o paired\_dat\_m180\_

#质控分析

fastqc -t 10 paired\_dat\_m3000\_1.fq paired\_dat\_m3000\_2.fq

fastqc -t 10 paired\_dat\_m180\_1.fq paired\_dat\_m180\_2.fq

#bowtie2比对两套高通量模拟的结果与参考基因组

bowtie2 -x ./Bowtie2Index/genome -U paired\_dat\_m180\_1.fq,paired\_dat\_m180\_1.fq -S ./test180.sam

bowtie2 -x ./Bowtie2Index/genome -U paired\_dat\_m3000\_1.fq,paired\_dat\_m3000\_2.fq -S ./test3000.sam

#samtools分析

#180

samtools stats ./data/test180.sam --ref-seq ./data/genome.fa \>data/sam.stats180.txt

samtools view -bS ./data/test180.sam \> ./data/test180.bam

samtools view ./data/test180.bam \> ./data/test2\_180.sam

samtools sort ./data/test180.bam -o ./data/test.sorted180.bam

samtools index ./data/test.sorted180.bam

samtools tview ./data/test.sorted180.bam

samtools depth ./data/test.sorted180.bam \>./data/bam.stat.depth180.txt

#3000

samtools stats ./data/test3000.sam --ref-seq ./data/genome.fa \>data/sam.stats3000.txt

samtools view -bS ./data/test3000.sam \> ./data/test3000.bam

samtools view ./data/test3000.bam \> ./data/test2\_3000.sam

samtools sort ./data/test3000.bam -o ./data/test.sorted3000.bam

samtools index ./data/test.sorted3000.bam

samtools tview ./data/test.sorted3000.bam

samtools depth ./data/test.sorted3000.bam \>./data/bam.stat.depth3000.txt

#测序深度数据挖掘

cut -f3 bam.stat.depth3000.txt |sort|uniq -c |sort -nr -k1 \> depth3000.txt

cut -f3 bam.stat.depth180.txt |sort|uniq -c |sort -nr -k1 \> depth180.txt

samtools flagstat /data/test180.bam

samtools flagstat /data/test3000.bam

#编写R语言脚本depth.r画depth180.txt和depth3000.txt的折线图

SOAPdenovo序列组装

配置文件lib.cfg

1
2
| max\_rd\_len=90[LIB]avg\_ins=180reverse\_seq=0asm\_flags=3rd\_len\_cutoff=90rank=1pair\_num\_cutoff=3map\_len=32q1=paired\_dat\_m180\_1.fqq2=paired\_dat\_m180\_2.fq[LIB]avg\_ins=3000reverse\_seq=1asm\_flags=3rank=2pair\_num\_cutoff=5map\_len=35q1=paired\_dat\_m3000\_1.fqq2=paired\_dat\_m3000\_2.fq |
| --- |

将程序挂上服务器，nohup.out作为记录工作日志

1
2
3
nohup SOAPdenovo-63mer all -s lib.cfg -K 29 -o SOAPdenovo\_out -p 20 &

makeblastdb -in GCA\_000975675.2\_Sc\_YJM1083\_v1\_genomic.fna -input\_type fasta -dbtype nucl -title YJM320 -out YJM320

1.3.3 同源搜索

从UniProt数据库下载所有已知真菌蛋白序列，关键词为taxonomy:fungi AND reviewed:yes，与基因序列数据fna文件进行tblastn比对，比对结果通过blast92gff3.pl程序生成注释数据gff文件。

运行代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#创建本地BLAST数据库

makeblastdb -in ./fugi\_review.fasta -input\_type fasta -title SC\_gDNA -dbtype nucl -out SC\_gDNA

#使用tblastn程序，把已知蛋白质序列和上述建立的本地BLAST数据库分别进行6和7比对

nohup tblastn -query ./fugi\_review.fasta -db SC\_gDNA -out ./blastx\_results.outfmt6 -evalue 1e-5 -outfmt 6 -max\_target\_seqs 1 -num\_threads 10 &

nohup tblastn -query ./fugi\_review.fasta -db SC\_gDNA -out ./blastx\_results.outfmt7 -evalue 1e-5 -outfmt 7 -max\_target\_seqs 1 -num\_threads 10 &

#运行blast92gff3.pl

perl blast92gff3.pl -exonType exon -geneType gene blastx\_results.outfmt6 \> pl\_output6\_.gff3

#利用gffcompare比对pl\_output6.gff3同基因组gff文件

./gffcompare -R -r GCA\_000975885.2\_Sc\_YJM320\_v1\_genomic.gff -o strtcmp pl\_output6.gff3

1.3.4 从头预测

利用augustus工具对基因序列数据fna文件进行从头预测, 提取输出文档中的有用信息进行同真菌蛋白序列进行blastp比对，从结果中获得具体的基因信息，构建基因注释数据gff文件。最后，利用gffcompare分别比对同源搜索结果、从头测序结果与原基因组

运行代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#创建本地BLAST数据库

makeblastdb -in ./fugi\_review.fasta -input\_type fasta -title uniprot\_protein -dbtype prot -out uniprot\_protein

#利用augustus对基因组进行从头预测

augustus --gff3=on --outfile=Sc\_augustus\_out.gff3 --species=saccharomyces\_cerevisiae\_S288C ./GCA\_000975885.2\_Sc\_YJM320\_v1\_genomic.fna

#进行blast比对

nohup blastp -query Sc\_augustus\_out.fa -out seq.blast -db uniprot\_protein -outfmt 6 -evalue 1e-5 -max\_target\_seqs 1 -num\_threads 20 &

#将比对结果seq.blast文件的预测基因名合并，使用R语言代码完成

#利用gffcompare将GCA\_000975885.2\_Sc\_YJM320\_v1\_genomic.gff与result.gff进行比对

gffcompare -R -r GCA\_000975885.2\_Sc\_YJM320\_v1\_genomic.gff -o strtcmp result.gff

#对实验三同源搜索所生成的进行同样的比对

gffcompare -R -r pl\_output6.gff3 -o pl result.gff

1.3.5 启动子分析与预测

从EPD数据库中下载任意一种启动子相关的DNA元件的HMM数据，去除序列中的小写atcg和未知N，对正反双链利用bootstrap法计算pvalue和score，取99%的置信区间绘制打分图和类ROC曲线（横坐标为位点，纵坐标为该距离内的启动子/所有启动子）

1.3.5 基因组可视化

使用IGV进行可视化。基因组序列：菜单栏 Genomes→Load Genome from File =》选中fasta 格式的基因组序列文档；同源预测（blast）结果：菜单栏 Files→Load from File=》选中blast 比对结果转换来的 gff3 文档；从头预测结果：菜单栏 Files→Load from File=》选中从头预测的gff3 文档;全部加载成功后，任选一个包含注释信息的区间截图保存。

二、结果

1. 基因组测序模拟

利用art illumine HighSeq2k对不同fold的数据分别模拟测序，结果见表2

表 2 测序数据与丢失值、覆盖率的关系

l f m 丢失值 覆盖率
100 1 150 3.70E-01 63.0%
100 2 150 1.37E-01 86.3%
100 3 150 5.05E-02 95.0%
100 5 200 7.42E-03 99.3%
80 10 150 4.76E-05 100% 
图 2 测序覆盖度f与覆盖率的关系

分析：在fold小于10的情况下，覆盖率随fold的增大而增大，其变化幅度随fold的增大而减小，f大于10后基本可保证完全覆盖。

2. 测序组装

2.1质控分析

分别对art illumina模拟180x和3000x的数据进行fastqc质控分析

图 3 fastqc 每个碱基位点的质量评估（√） 图 4 fastqc平均序列质量（分布近似正态√）

图 5 fastqc GC含量（理论与实际相近） 图 6 fastqc 序列重复度（√）

其峰值对应GC含量与基因组相近（约38）

Fastqc结果数据符合质控标准，通过质控检验。

2.2数据比对

利用bowtie2分别进行两套数据同参考基因组的比对和samtools分析，这里我特别分析了测序深度方面的数据。结果如下

图 7 flagstat 180x

图 8 flagstat 3000x

二者的测序均在98%以上，情况符合预测。

图 9 测序深度与数量关系 上面为180,下面为3000

趋势正常（测序深度集中在45左右，符合45x），大致符合正态分布。

2.3序列组装

利用SOAPdenovo软件对序列片段进行序列组装具体结果汇总表如下：

Table 1 序列组装结果数据

数量 N50
Scaffold 265 258107
contig 43082（大于100bp） 273

2.4覆盖率计算

利用网络平台QUAST进行可视化及实际覆盖率计算。

图 10 Quast可视化

图 11 Quast数据

利用quast分析所得的覆盖率为87.6%。

3 同源搜索

将已知真菌蛋白序列与基因序列数据fna文件进行tblastn比对，比对结果通过blast92gff3.pl生成注释数据gff文件。

4. 从头预测

利用augustus工具针对酵母基因组序列数据进行从头预测，并将结果文档同真菌蛋白序列进行blastp比对，从结果中获得具体的基因信息，补全基因注释数据gff文件。

利用gffcompare分别比对同源搜索结果、从头测序结果与原基因组,通过总结画出同源搜索、从头预测、原始基因组的对比表

表 3 同源搜索、从头预测、原始基因组的gffcompare

从头预测的结果外显子、内含子等各项情况均远优于同源建模，有理由认为从头预测优于同源建模。

5 启动子分析与预测

根据TATA box的HMM数据，对正反双链利用bootstrap法计算pvalue和score，分别绘制1号和12号染色体的打分图

图 12 染色体1全部序列的打分图

图 13 染色体12全部序列的打分图

从12号染色体的打分图中，我们也可以看出，该染色体有很长一段区域打分情况一直稳定在一个固定的值处，可以认为并没有被匹配上TATA box序列

绘制类ROC曲线（横坐标为位点，纵坐标为该距离内的启动子/所有启动子）

图 14 染色体1 +链TATAbox

图 15 染色体1 -链TATAbox

匹配得到的TATA box序列大多在基因前5000bp内

6. 基因组可视化

图 16 基因组可视化结果

从上到下依次为原基因组、同源建模结果、从头预测结果、TATA box预测结果可视化结果。